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ABSTRACT 

We introduce two new methods to obtain reliable velocity field statistics from N-body simula- 
tions, or indeed from any general density and velocity fluctuation field sampled by discrete points. 
These methods, the Voronoi tessellation method and Delaunay tessellation method, are based on 
the use of the Voronoi and Delaunay tessellations of the point distribution defined by the locations 
at which the velocity field is sampled. In the Voronoi method the velocity is supposed to be uniform 
within the Voronoi polyhedra, whereas the Delaunay method constructs a velocity field by linear 
interpolation between the four velocities at the locations defining each Delaunay tetrahedron. 

The most important advantage of these methods is that they provide an optimal estimator 
for determining the statistics of volume-averaged quantities, as opposed to the available numerical 
methods that mainly concern mass- aver aged quantities. As the major share of the related analytical 
work on velocity field statistics has focussed on volume-averaged quantities, the availability of 
appropriate numerical estimators is of crucial importance for checking the validity of the analytical 
perturbation calculations. In addition, it allows us to study the statistics of the velocity field in 
the highly nonlinear clustering regime. 

Specifically we describe in this paper how to estimate, in both the Voronoi and the Delaunay 
method, the value of the volume-averaged expansion scalar 9 = H^V ■ v - the divergence of the 
peculiar velocity expressed in units of the Hubble constant H, as well as the value of the shear and 
the vorticity of the velocity field, at an arbitrary position. The evaluation of these quantities on 
a regular grid leads to an optimal estimator for determining the probability distribution function 
(PDF) of the volume-averaged expansion scalar, shear and vorticity. Although in most cases both 
the Voronoi and the Delaunay method lead to equally good velocity field estimates, the Delaunay 
method may be slightly preferable. In particular it performs considerably better at small radii. 
Note that it is more CPU-time intensive while its requirement for memory space is almost a factor 
8 lower than the Voronoi method. 

As a test we here apply our estimator on an N-body simulation of such structure formation 
scenarios. The PDFs determined from the simulations agree very well with the analytical predic- 
tions. An important benefit of the present method is that, unlike previous methods, it is capable 
of probing accurately the velocity field statistics in regions of very low density, which in N-body 
simulations are typically sparsely sampled. 

In a forthcoming paper we will apply the newly developed tool to a plethora of structure 
formation scenarios, both of Gaussian and non-Gaussian initial conditions, in order to see in how 
far the velocity field PDFs are sensitive discriminators, highlighting fundamental physical differences 
between the scenarios. 

Subject headings: Cosmology: theory - large-scale structure of the Universe - Methods: numerical 
- statistical 

1. Introduction 

Besides the distribution of galaxies and the temperature fluctuations in the cosmic microwave 
background, and possibly weak gravitational lensing of background galaxies by large scale structures 
(see e.g. Villumsen 1995), the velocity field on cosmological scales is one of the main sources of 



information on the formation and evolution of structure in the Universe. Early work indicated 
the existence of large-scale velocity flows (Rubin et al. 1976) and established the existence of the 
motion of the Local Group with respect to the rest-frame defined by the microwave background 
(Smoot & Lubin 1979). However, it was the work by Burstein et al. (1987) that established beyond 
doubt that the Local Group is participating in a large scale streaming motion. 

The advent of reliable redshift-independent distance estimators lead to an enormous growth of 
activity in the field of measuring and interpreting the peculiar velocities of galaxies. This growth 
of attention was evidently fed by the fact that the velocity field provides direct information on the 
dynamics of the Universe on scales of more than a few Mpc. Above these scales dynamical relaxation 
has not had yet a chance to wash out the memory of the conditions in the early Universe. The 
velocity field can therefore be fruitfully investigated by means of perturbation analysis. Particularly 
useful is the, O-dependent, velocity-density relationship that follows from linear theory (see e.g. 
Peebles 1980). Moreover, a general result of perturbation theory predicts that the rotational part 
of the velocity field vanishes. Based on this observation Bertschinger & Dekel (1989) developed the 
non-parametric POTENT method in which the local cosmological velocity field is reconstructed 
from the measured line-of-sight velocities. In a series of papers (e.g. Bertschinger et al. 1990), they 
applied their method to the existing catalogues of galaxy peculiar velocities. 

The POTENT analysis also paves the way for a more quantitative analysis of the velocity field, 
estimating various statistical properties and their relation to the properties of the density field. 
For instance, via the velocity-density relationship it is possible to reconstruct the corresponding 
density field or, by comparison with a uniformly sampled galaxy redshift catalogue (in particular 
the IRAS based redshift catalogues, e.g. Strauss et al. 1990), obtain an idea of in how far the 
galaxy distribution forms a biased tracer of the underlying mass distribution. If it is assumed that 
this bias can be simply represented by a linear bias factor b, the value of fi a6 /6 can be estimated 
from the measured peculiar velocity field and a uniformly sampled redshift catalogue (e.g. Yahil et 
al. 1991). 

A variety of other methods have been proposed to determine this combination of and b (see 
the review paper of Dekel 1994 and references therein). On the basis of a more specific statistical 
analysis, in particular of the velocity divergence, other studies managed to determine these two 
parameters separately. For example, Nusser & Dekel (1993) proposed to use a reconstruction 
method assuming Gaussian initial conditions to constrain Q, while Dekel & Rees (1994) used voids 
to achieve the same goal. 

Another approach has been proposed by Bernardeau et al. (1995) and Bernardeau (1994a). 
This is based on the use of statistical properties of the divergence of the locally smoothed velocity 
field, its skewness (the 3 rd order moment) and its kurtosis (4 th order moment). It was shown that 
these statistical quantities are potentially very valuable to measure f2 or to test the gravitational in- 
stability scenario. The early analytical work was subsequently extended towards the determination 
of the complete velocity divergence probability distribution function for gravitational instability 
scenarios starting from Gaussian initial conditions in the case the velocity field is filtered by a top- 
hat window function. Preliminary comparisons with numerical simulations (Bernardeau 1994b, 
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Juszkiewicz et al. 1995, Lokas et al. 1995) yielded encouraging results. However, a comparison of 
the analytical results with the PDFs determined from N-body simulations is complicated consider- 
ably by the fact that in the latter velocities are only known at, non-uniformly distributed, particle 
locations. 

In this paper we address specifically the issue of the discrete nature of the velocity sampling, 
which leads to the development of a numerical method to obtain reliable velocity field statistics 
from N-body simulations. The goal is threefold. First of all, we wish to have an independent 
way of checking whether the perturbation calculations that yield the quasi-linear results are indeed 
valid. Secondly, if so, these analytical results form a good 'calibration' point for the numerical tool 
that we have developed here, so that we can use it with confidence in highly nonlinear conditions. 
Finally, because the velocity field in observational samples is also only known at a discrete number 
of positions, the location of galaxies, we may be able to apply the developed method, in adapted 
form, to the available catalogues of measured galaxy peculiar velocities. 

The central problem that we address here is that while almost all analytical results concern 
volume- averaged quantities, almost all available numerical estimators in essence only yield mass- 
averaged quantities. This may considerably complicate any comparison, and even lead to false 
conclusions regarding e.g. the validity of perturbation theory. To improve upon this situation we 
introduce two new numerical methods, the Voronoi tessellation method and the Delaunay tessella- 
tion method. Both methods are based on two important objects in stochastic geometry, the Voronoi 
and the Delaunay tessellation of the point set consisting of the points at which the velocity field has 
been sampled. A Voronoi tessellation of a set of nuclei is a space-filling network of polyhedral cells, 
each of which delimit that part of space containing that is closer to its nucleus than to any of the 
other nuclei. The Delaunay tessellation is also a space-filling network of mutual disjunct objects, 
tetrahedra in 3-D. The four vertices of each Delaunay tetrahedron are nuclei from the point set, 
such that the corresponding circumscribing sphere does not have any of the other nuclei inside. 
The Voronoi and the Delaunay tessellation are closely related, and are dual in the sense that one 
can be obtained from the other. 

The earliest use of Voronoi tessellations can be found in the work of Dirichlet (1850) and 
Voronoi (1908) in their work on the reducibility of positive definite quadratic forms. However, their 
application to random point patterns has caused this concept to arise independently in various fields 
of science and technology, ranging from molecular physics to forestry (see Stoyan, Kendall, and 
Mecke 1987 for references). Because of these diverse applications they acquired a set of alternative 
names, such as Dirichlet regions, Wigner-Seitz cells, and Thiessen figures. Despite the simplicity 
of its definition, analytical work on the statistics of Voronoi tessellations has appeared to be rather 
complicated and cumbersome, so that present analytical knowledge is mainly restricted to a few 
statistical moments of the distribution function of geometrical properties of Voronoi tessellations 
generated by homogeneous Poisson processes (see e.g. Meyering 1953, Gilbert 1962, Miles 1970, 
and M0ller 1989). For the time being, numerical work is therefore an inescapable necessity for any 
progress in this field (see e.g. Van de Weygaert 1991b, 1994). 

Within astronomy, and in particular cosmology, Voronoi tessellations are mostly known for 
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their application as geometrical models for astrophysical structures. One of the first applications 
was by Kiang (1966), who used them in an attempt to derive a mass spectrum resulting from the 
fragmentation of interstellar clouds. The similarity of Voronoi tessellations to the cellular patterns 
in the galaxy distribution sparked a lot of work in cosmology. Matsuda & Shima (1984) pointed 
out the similarity between two-dimensional Voronoi tessellations and the outcome of numerical 
clustering simulations in a neutrino-dominated universe. Independently, Voronoi tessellations were 
introduced into cosmology in a study by Icke & Van de Weygaert (1987) of the statistical properties 
of two-dimensional Voronoi tessellations. They argued that a cellular pattern similar to Voronoi 
tessellations is a natural outcome of the evolution of a void-dominated Universe. Their statisti- 
cal analysis was extended to three dimensions in Van de Weygaert (1991b, 1994), based on the 
completion of a three-dimensional geometrical Voronoi algorithm. In an early analysis of three- 
dimensional Voronoi tessellations, Van de Weygaert & Icke (1989) found that Voronoi tessellations 
also possess some interesting clustering properties, as was confirmed in a Monte-Carlo study by 
Yoshioka & Ikeuchi (1989). Since then, the amount of applications of Voronoi tessellations as a 
useful, conceptually simple model for a cellular or foam- like distribution of galaxies on large scales, 
has steadily increased (Coles 1990; Van de Weygaert 1991a,b; Ikeuchi & Turner 1991; SubbaRao & 
Szalay 1992; Williams 1992; Williams et al. 1992; Goldwirth, Da Costa, k Van de Weygaert 1995). 

However, their use as a geometrical model is a different class of applications of Voronoi tes- 
sellations than the one we use in the present paper. Here we follow another philosophy, namely 
that the sensitivity of the geometrical characteristics of the Voronoi tessellations to the underly- 
ing nucleus distribution makes the Voronoi tessellation, and its dual the Delaunay tessellation, a 
potentially very useful instrument to study the properties of a point process. Their usefulness 
as statistical descriptors was suggested earlier by e.g. Finney (1979), who introduced the name 
'polyhedral statistics'. Instead of being interested in the generating point process itself, we turn 
our attention to developing a reliable and/or optimal description of the velocity field sampled by 
the point process. 

The first method that we have introduced here is based on the Voronoi tessellations, and 
follows directly from the definition of volume-averaged velocities (eqn. 3). It can be considered as 
the multidimensional extension of the approximation of a function of one variable by a constant 
value in a finite number of bins, the constant value being equal to the function value at the point in 
the bin. The Voronoi method yields a velocity field with constant values of the velocity components 
within each Voronoi cell of the defining point distribution. The velocity within the whole interior of 
the cell is equal to the velocity of the nucleus of that cell. Consequently, only at the boundaries of 
the Voronoi cells the velocity gradient has a non-zero value. The subsequent operation of volume- 
averaging of a quantity therefore consists of determining the intersection of Voronoi walls with the 
appropriate filter, for a top-hat filter a sphere of radius R, and weighing the value of that quantity 
in the wall with the size of intersection area. 

The Delaunay method, on the other hand, should be regarded as the multidimensional recipe 
for linear interpolation. In a space of dimension d > 2 linear interpolation consists of assuming 
constant function gradients in interpolation intervals defined by d + 1 points, 'hyper-triangles'. 



4 



In two-dimensional space a hyper-triangle is a triangle, in three-dimensional space a tetrahedron. 
Unlike the one-dimensional case, the choice of multidimensional interpolation tetrahedra may not 
be unique. However, here we argue that Delaunay tetrahedra are a natural and logical choice 
based on the requirements that the whole of the sample space is covered by a space-filling network 
of mutually disjunct tetrahedra and that these tetrahedra should be as compact as possible to 
minimize approximation errors. The (constant) value of each of the 9 velocity field gradients in 
each of the Delaunay tetrahedra is determined from the locations of and velocities at the four points 
that define each of them. Summarizing, the Delaunay method consists of three major steps: (1) 
construction of the Delaunay tessellation, (2) determination of the 9 velocity gradients dvi/dxj 
in each of the tetrahedra, (3) volume averaging of the obtained field of velocity gradients. For a 
top-hat filter the latter step consists of determining the volume of the intersection of the Delaunay 
tetrahedra with the filter sphere. 

In this paper, we start by shortly discussing the conventional methods to sample the value of 
velocity fields on grid positions from the value of the velocity at a discrete number of points. In 
particular we stress the fundamental difference between mass and volume weighted velocity av- 
erages. Conventional estimators are almost always based on mass weighted velocity fields. This 
may induce complications in the comparison with theoretical predictions. This leads to the in- 
troduction in section 3 of the Voronoi and the Delaunay method. Both are good estimators for 
volume weighted velocity fields, and are therefore instrumental in improving comparisons between 
theoretical predictions and the results of N-body simulations. The accompanying appendices A 
and B contain detailed descriptions of the geometrical details of these sampling techniques. Com- 
putational considerations are discussed in section 4, while section 5 contains a discussion of an 
application of both methods to an N-body simulation of galaxy clustering, determining the values 
of the divergence, vorticity and shear of the velocity field on a regular grid. These values are com- 
pared with each other as well as with the ones determined by a conventional estimator. Finally, 
in section 6 we conclude with a discussion of the virtues of the new methods, and suggestions for 
future applications. 

2. Discretely sampled velocity field 

The fact that the velocity field is only known at a finite number of discrete positions is a major 
technical obstacle for obtaining reliable estimates of statistical parameters of the velocity field. It 
is of importance both in the observational data as well as in N-body simulations. One possibility 
is to smooth the galaxy velocity field with a filter. For example, Bertschinger et al. (1990) choose 
to filter the measured galaxy velocities with a Gaussian smoothing function. In one case they took 
a fixed smoothing length, to be preferred for a rigorous statistical analysis, while for an optimal 
representation of the velocity and density field they adopted a filter with an adaptive smoothing 
length. By taking this length to be equal to the distance to the fourth nearest neighbour the 
analysis automatically becomes better in the well-sampled high-density regions, thereby reducing 
the problem of noisy data and sparsely sampled underdense regions. In their analysis of numerical 
simulations Juszkiewicz et al. (1995) and Lokas et al. (1995) also used smoothing of the velocity 
field by a Gaussian filter with a fixed smoothing length to obtain the local velocity field on a regular 
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grid. 

Effectively these are mass weighted velocity fields, 

J dx v(x) p(x)Wm(x, x ) 

V m ass(xo) = J , (1) 

where Wm(x, xo) is the used filter function defining the weight of a mass element in a way dependent 
on its position relative with respect to the position xo. In other words, v mass is effectively the 
velocity corresponding to the average momentum within the filter volume. For a discrete particle 
distribution v mass reduces to, 

J dxv(x) Wm(x,x ) ^5d(x - Xj) 
/ dx^ M (x,xo)^5 jD (x-x i ) 

(2) 

^WjV(Xj) 



E 



where it>j = WM(xj,xo) and <5d(x, xo) the Dirac delta function. A major complication of such 
an analysis is that a comparison of statistical properties of the velocity field obtained in this way 
with the known analytical results is not straightforward. Almost without exception these analytical 
results are based on the volume weighted filtered velocity field v, 

/ dxv(x)Wy(x,xo) 
v(x ) = J— r , (3) 

dx Wy(x, Xq) 



where VFy(x, xo) is a possibly used weight function. At the present the only analytical work that 
has treated certain aspects of the statistics of the mass averaged velocity, namely the skewness 
of 6 = H~ l V ■ v, is the work by Bernardeau et al. (1995). A major complication is that this 
quantity involves the density field, and therefore would introduce the unknown bias between mass 
and galaxy density field in a practical implementation. 

In order to get an estimate of the volume averaged velocity Juszkiewicz et al. (1995) use a two- 
step scheme, wherein they first determine a velocity on a grid according to equation (2), and then 
use the resulting grid of velocities to determine volume averaged velocities according to equation 
(3), the volume averaging being accomplished by Gaussian filtering. Bernardeau (1994b) followed a 
similar procedure, whereby he used a top-hat filter for the volume averaging. Such a scheme would 
yield reliable results if the filter length of Wm would be much smaller than that of Wy . However, 
for technical reasons this is often difficult to attain. For example, it requires a very small grid size 
which would be excessively computer time and memory consuming. 

In the study of Lokas et al. (1995) the same approach was followed in a comparison between 
N-body simulations and analytical perturbation calculations. While they found good agreement 
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Figure 1. The Dclaunay triangulation (left frame) and Voronoi tessellation (right frame) of a distribution of 25 nuclei (stars) 
in a square (central panel). Periodic boundary conditions are assumed. 

with the obtained density field moments, there already appeared to be substantial discrepancies 
in the case of the velocity divergence 9 by the time oq ~ 0.1. They suggested as possible causes 
for the disagreement that (1) the perturbation approximation already starts to break down by the 
time as i=3 0.1, (2) the N-body simulations do not evolve the velocity divergence field with sufficient 
accuracy and (3) the estimator of the smoothed velocity divergence from the N-body simulations is 
too noisy or too inaccurate. In order to improve upon this situation it is therefore highly desirable 
to develop a more reliable estimator, preferentially closer related to the 'volume-averaged' nature 
of the perturbation calculations. In the following sections we will introduce two new methods that 
overcome the dilemma of a required very small initial smoothing length by constructing the Voronoi 
and Delaunay tessellations of the point distribution. 

3. Voronoi and Delaunay Tessellations 

The first new estimator that we introduce here is based on Voronoi tessellations. It follows in 
a rather direct way from an asymptotic interpretation of the definition of mass filtered quantities 
(eqn. 2). Although it leads to good results, it corresponds to an artificial situation of a discon- 
tinuous velocity field. This is successfully improved upon by a subsequent further elaboration and 
extension of the Voronoi method to a second estimator, based on the division of space into Delaunay 
tetrahedrons. 

3.1 The Voronoi Tessellation 

In section 2 we made the observation that a good approximation of volume averaged quantities is 
obtained by volume averaging over quantities that were mass filtered with, in comparison, a very 
small scale for the mass weighting filter function. We can then make the observation that the 
asymptotic limit of this, namely using a filter with an infinitely small filter length (see eqn. 2), 



*( x o) 



i 

N 



Wi 

i=2 



v(xx) +]T^-v( Xi ) 



N 

Wi 



(4) 

v(xi) , 



where we have ordered the locations i by increasing distance to xo and thus by decreasing value of 
Wi, is in fact taking the velocity v(xi) of the closest particle as the estimate of the mass averaged 
velocity. In other words, we can divide up space into regions consisting of that part of space closer 
to a particle than to any other particle, and taking the velocity of that particle as the value of the 
velocity field in that region. The search for the closest particle of each point of the field naturally 
leads to the construction of the Voronoi tessellation associated with the particle distribution. This 
division of space is a familiar and important concept in the field of stochastic geometry (see Stoyan, 
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Figure 2. Three examples of three-dimensional Voronoi cells. They are taken from a Voronoi tessellation generated by 1000 
Poissonian distributed nuclei in a box with periodic boundary conditions. The stars in each of the cells represent the position 
of the cell nucleus. 

Kendall and Mecke 1987 for an overview of this field). It consists of space-covering and mutual 
disjunct set of convex cells, each of which contains a particle of the original particle distribution, 
enclosing the points of space for which the closest particle is precisely the one in the cell. 

Formally a Voronoi tessellation (Voronoi 1908, Dirichlet 1850) can be defined as follows. As- 
sume that we have a distribution of a countable set <& of nuclei {xi} in Let X\,X2,Xs, ... be the 
coordinates of the nuclei. Then, the Voronoi region Hj of nucleus i is defined by the following set 
of points x of the space: 

LTj = {x\d(x,Xi) < d(x,Xj) for all j ^ i}, (5) 

where d(x, y) is the Euclidian distance between x and y. In other words, IT is the set of points which 
is nearer to Xj than to Xj, j ^ i. Each region IL, therefore consists of the intersection of the open 
half-spaces bounded by the perpendicular bisectors of the segments joining x\ with each of the other 
Xj's. Hence, Voronoi regions are convex polyhedra (3-D) with finite size according to definition (5). 
Each Ilj is called a Voronoi •polyhedron. The complete set of {II j} constitute a tessellation of 
the Voronoi tessellation V(<1?) relative to <I>. A two-dimensional Voronoi tessellation of 25 cells is 
shown in figure 1 (right-hand frame), while an idea of a three-dimensional Voronoi tessellation can 
be obtained from the three Voronoi cells displayed in figure 2. The latter three cells are taken from 
a Voronoi network of 1000 cells, generated by Poissonian distributed nuclei, in a box with periodic 
boundary conditions. 

As described above, given a field of velocities at a set of discrete particles a reasonable first 
assumption is that the velocity is constant within each Voronoi cell. The Voronoi method can there- 
fore be regarded as the three- (or two-) dimensional equivalent of approximating a one-dimensional 
function f(x) by a sequence of intervals wherein the function has a constant value. If the value 
of the function f(x) is known at M discrete points Xi (upper panel figure 4), this one-dimensional 
equivalent of the Voronoi method consists of adopting a constant value f(xi) in the interval between 
(xi + Xj_i)/2 and (xi + Xj+i)/2 (see central panel of figure 4). 

The first step of the Voronoi algorithm consists of calculating the Voronoi tessellation that 
is defined by the set of points at which the velocity field has been sampled. For this we use the 
three-dimensional geometrical Voronoi code that was developed by Van de Weygaert (1991b, 1994). 
Starting from an input of points this code calculates the complete geometrical structure, i.e. the 
location of the walls, edges and vertices, of the corresponding Voronoi tessellation. 

Subsequently, we proceed by determining the corresponding volume averaged quantities. This 
is accomplished by a volume filtering of the resulting velocity field. By adopting a top-hat filter 
Wth w hh radius R as the volume filter Wy , the problem of determining the corresponding volume- 
averaged velocity gradient Vij = VjVi is determining the average value of Vij within a sphere of 
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Figure 3. Two-dimensional illustration of the Voronoi method to estimate velocity field gradients (compare eqn. 7-9). The 
stars indicate the positions of the nuclei of the Voronoi cells. Part of the cell boundaries are indicated by the walls with finite 
thickness. The dashed arrows at each of the stars are the velocity vectors at each of these locations. The solid arrow is the 
normal vector of the wall k, of thickness Asj., between the two cells fci and k2, in whose interior the velocity field has values 
Vki and Vk2 respectively. Partially indicated is the circle corresponding to a top-hat filter centered on a grid position. 

radius R, 

dXj J dx^ Ti/ (x,x ) (6) 

3 f A 

fllX Vj 



^j( x o) = "cTT" = 



J R ^ ' 

where the latter volume integral is over the part of space enclosed by the sphere with radius R 
centered on xo. The constant value of the velocity v within each Voronoi cell automatically implies 
that the value of the velocity gradient fjj(x) is equal to zero in their interior, so that the cells 
themselves have a contribution zero to the integral / dxvij. Only at the boundaries between the 
Voronoi cells Vij will have a non-zero value. Moreover, the value of Vij will be constant within each 
wall, as it corresponds to the change of value of v between the two corresponding neighbour cells 
(see figure 3). The finite contribution by any Voronoi wall that lies within or intersects the filter 
sphere can be calculated by considering the volume defined by the surface of the wall and having 
an infinitesimal width As perpendicular to the wall. Imagine a Voronoi wall k, being the boundary 
between two Voronoi cells kl and k2 in whose interior the velocities are v&i and Vfc 2 (see figure 
3), that intersects the top-hat sphere. The wall has a perpendicular width As k and a surface area 
Ak within the sphere, while its orientation is determined by its normal vector n^. The velocity 
gradient v XyV can then be easily inferred from the velocity change Av x = (v& 2 — Vfci) • e x along the 
x-direction. Because this corresponds to an interval Ay = (n^ • e y ) As^ in the y direction, we find 

dv x _ A^£ = (n fc • e y ) (v fc2 - Vfci) • e x 
dy ~ Ay ' As k 

This can be generalized directly to yield the following result for the volume averaged value v^j 
(with = 1,2,3), 

5 *J = T~BS 12 A k ( n k ■ ej) (v fc i - v fc2 ) • ej , (8) 

where the sum is over all walls that intersect the sphere. In the specific case of the volume-averaged 
velocity divergence 6, this leads to the expression 

~ V • v 3 ^ A k n k ■ (vfci - v fc2 ) /n , 

9 = ^r = ^¥\ h • (9) 

The problem of determining the volume-averaged velocity velocity gradients (eqn. 8), or velocity 
divergence (eqn. 9), has therefore been reduced to determining the intersection of the walls in the 
Voronoi tessellations with spheres. This is a geometrical problem that can be solved with some 
effort (see Appendix A). 

In practice we repeat this procedure of top-hat averaging at each point of a regular grid. 
From the values of dvi/dxj at each grid point we can then easily evaluate the value of the velocity 
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Figure 4. Illustration of the one-dimensional equivalents of the Voronoi and the Delaunay method. The value of some function 
f(x) has been determined at 20 randomly distributed positions Xi along a line (top panel). The height of a star at position 
Xi is proportional to the value f(xi). Central panel: 'Voronoi method'. The function / in the whole bin centered on a point 
Xi, bounded by the points halfway in between X{ and its neighbours X{—i and X{+ i, is assumed to have the constant value 
f(xi). Note that these bins are essentially one-dimensional Voronoi cells. Bottom panel: 'Delaunay method'. The function 
/ is approximated by linear interpolation from f(xi) to /(aii+i) in the interval (xi, Zi+i), and similarly between f(xi) and 
f(xi-i) in the interval (xi—i,Xi). The intervals (xi,Xi+i) and (xi-i,Xi) can be considered as the one-dimensional equivalents 
of Delaunay cells. 

divergence 6, the shear <7jj and the vorticity u, where w = Vxv = e^Wij (and e kl ^ is the completely 
antisymmetric tensor), 



Summarizing, the end result of the operations described above consists of fields of top-hat averaged 
quantities like velocity divergence, shear and vorticity, sampled at regular grid intervals. 

3.2 The Delaunay Tessellation 

A major characteristic of the Voronoi approach is that it leads to a discontinuous velocity 
field. This is evidently not the only unique way of defining the velocity field from a discrete set 
of sample points. Moreover, the Voronoi method may have problems in the case of small filter 
radii. Relevant non-zero values for the velocity gradients are only produced in the Voronoi walls. 
Many filter spheres would remain empty when the filter scale is smaller than the average distance 
between Voronoi walls, i.e. when their scales are smaller than ps L/N 1 ^ 3 (with L the boxsize and 
N the number of Voronoi cells). This would yield many irrelevant and noisy filter velocity gradient 
averages. It may therefore be worthwhile to define another independent method based on a different 
interpolation scheme. In fact, this would give us the possibility of internally checking the results of 
the new methods that we have introduced here. 

Indeed the Voronoi method can be viewed as an elementary zeroth-order interpolation scheme. 
However, it is possible to define a velocity field based on linear interpolation between the velocities 
of the sample points. This is the multidimensional equivalent of the one-dimensional situation where 
the approximation of a function by constant function values in bins centered on the sample points 
(central panel figure 4) is replaced by linearly interpolated function values in between those sample 
points (lower panel figure 4). For d = 1 linear interpolation simply consists of the approximation 
of a function f(r) by the value f(r) = a.i f(ri) + CKj+i /(rj+i) in the interval V{ < r < rj+i, where 
< a < 1 and ccj + CKj+i = 1. The natural extension to a space of arbitrary dimension d of the 
concept of the one-dimensional interpolation interval ri < r < rj+i is a 'hyper-triangle' defined by 
d+1 vectors I"*,, the vertices of that object. For d = 2 this is a triangle, for d = 3 a tetrahedron. Any 
vector r within the 'hyper-triangle' is a linear combination of the d+1 vectors r^, r = J2t=i a k r k > 
with < ak < 1 and a± + . . . + a^+i = 1. The linearity of the approximation of the function /(r) 
then implies that 




(10) 




dvi dvj 
dxj dxi 
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/(r) = /( ri ) + a 2 V/-(r 2 - ri ) + ... + 

+ ad+i V/ • (r d+1 - ri) 
= /(r 1 )+« 2 (/(r 2 )-/(r 1 )) + ...+ 

+ a d+ i (/(r d+ i) - /(ri)) 

d+1 
fc=l 

The choice of appropriate, and if possible optimal, interpolation 'hyper-triangles' is a critical issue 
that to a considerable extent determines the quality and accuracy of the multidimensional linear 
interpolation. An obvious minimal requirement is that the linear approximation of the function / 
leads to a unique value at every point in the subspace defined by the set V . In other words, we 
need a space-filling covering by mutual disjunct tetrahedra (for now we will restrict ourselves to 
d = 3, although it is trivial to follow the same argument for any other d). In the interior of each 
of these tetrahedra each of the velocity gradients dvi/dxj has one particular constant value, whose 
value is determined by the value of the function / at the locations of its four vertices. A crucial 
requirement for an optimal accuracy of the linear approximation (eqn. 11) is that the tetrahedra 
in the space-filling network are as compact as possible, in the sense of having a size and elongation 
that are as small as possible. A uniquely defined solution for such an optimal 'triangulation' does 
not exist or, rather, is not really known. Here we argue that a Delaunay tessellation (Delone 1934) 
is at least a good approximation of such an optimal triangulation. 

Formally, the Delaunay tessellation V(P) of a point set V is defined to be the tessellation 
consisting of all the tetrahedra defined by four nuclei whose circumscribing sphere is empty in 
the sense that no nucleus of the generating set of nuclei should be inside the circumsphere. A 
two-dimensional illustration of a Delaunay triangulation is given in figure 1 (left-hand frame). A 
principal characteristic of the circumsphere of a Delaunay tetrahedron is that the centre of the 
circumsphere (circumcentre) is a vertex of the Voronoi tessellation. This follows immediately from 
the extrapolation of the definition. After all, according to the definition of the Voronoi tessellation 
a vertex is defined by four nuclei that are equidistant from the vertex, i.e. the vertex is the 
circumcentre of the circumsphere of these four nuclei. If then there were a fifth nucleus within the 
sphere, it would be nearer to the circumcentre than the four nuclei on the surface of the sphere. 
This would imply that the centre cannot be the common vertex of the Voronoi polyhedra of these 
four nuclei. Ergo, the four nuclei have to define a Delaunay tetrahedron. In other words, the 
Delaunay tessellation is the network that is obtained by joining all pairs of nuclei in V whose 
Voronoi polyhedra share a Voronoi wall. Such a pair of nuclei is called a contiguous pair. The close 
relationship between the Delaunay and the Voronoi tessellation can be quite well appreciated from 
the right-hand frame of Figure 1, depicting the Voronoi tessellation of the same point set as the 
one in the left-hand frame. 

Besides the fact that Delaunay tetrahedra fulfil the requirement of compactness, them being 
objects of minimal size and elongation, an additional and equally important argument to the 
use of the Delaunay tetrahedra as the choice for linear interpolation intervals is provided by the 



11 



Figure 5. A visual impression of the Voronoi (left-hand panel) and the Dclaunay method (right-hand panel) for approximating 
a function /(r) = f{r\,T2) of two variables r\ and r%- Left-hand panel: the Voronoi method. The function / has been measured 
at the location of the black dots. In the (ri,T2) plane the corresponding Voronoi cells (neighbours of each other) have been 
indicated by solid lines. The resulting function field is indicated by the lightly shaded polygons, the Voronoi cells around the 
central nuclei r;. Each of these polygons is positioned at the corresponding height f(ri). For comparison with the Dclaunay 
method, the darkly shaded triangle indicates the value of the function field inside one of the corresponding Dclaunay triangles. 
Right-hand panel: the Dclaunay method. The dots indicate the positions r; of the points at which the function /(ri) has been 
determined. The same five Voronoi cells as in the top panel are indicated by the solid lines in the (r\,T2) plane, while all the 
Dclaunay triangles corresponding to these cells have been indicated by the dashed lines. The resulting function field is one of 
differently inclined triangles (shaded), connecting with each other at the boundaries of the Delaunay triangles. Evidently, at 
each of the locations ri (black dots) the function / has the value 



duality between Delaunay and Voronoi tessellations. Linear interpolation requires the definition of 
neighbour intervals. An arguably natural definition of 'neighbour points' in the multidimensional 
situation is that the two points share a Voronoi wall, i.e. they should be contiguous to each other. 
Delaunay tetrahedra therefore seem to be a natural choice for linear interpolation intervals defined 
by four nuclei. 

Moreover, by virtue of their duality, the Voronoi and the Delaunay tessellation for a given point 
set V are calculated by the same three-dimensional geometrical Voronoi tessellation code that was 
mentioned in the previous subsection (Van de Weygaert 1991b, 1994). In this case the output is 
limited to a listing of all vertices of the Voronoi tessellations and the location of the four generating 
points that define the corresponding Delaunay tetrahedron. 

To give a visual impression of the relationship between the Voronoi and the Delaunay approx- 
imation method, figure 4 shows the resulting values of a function /(r) = jiv^r-i) for d = 2. The 
Voronoi method (left-hand frame) yields a field that consists of regions of a constant value, Voronoi 
cells in (r±,r2) plane. The resulting image is therefore one of a field of pillars of different height, 
with a Voronoi cell at the base of each of the pillars. The Delaunay method divides space into 
triangular regions, the Delaunay triangles in the (ri, r2) plane, in which not the field value but the 
field gradients are constant. This yields a field of differently oriented triangles, connecting at their 
edges to the neighbouring triangles (right-hand frame figure 5). 

Having defined the Delaunay tetrahedra as the interpolation intervals, we proceed by deter- 
mining the values of the (constant) values of each of the nine velocity gradient tensor components 
dvi/dxj in each of these tetrahedra. These are determined from the location of each of the four 
vertices, ro, ri, T2 and r3, and the value of the velocity field at each of these locations, vo, Vi, V2 
and V3. Defining the quantities Ax n = x n — xo, Ay n = y n — Vo and Az n = z n — zo, for n = 1, 2, 3, 
as well as Av xn = v xn — v x o, Av yn = v yn — v y o and Av zn = v zn — v z o, the following nine linear 
relations are obtained, with n = 1,2,3, 



dv x dv x dv x 

Av xn = -T—Ax n + — — Ay n + — — Az n 
ox oy oz 



Av 



dv 



d 



t/n - -rr^-Axn + —^-Ay n + — —Az n 

ox oy oz 

dv z dv z dv z . 

Av zn = -r-Ax„ + -^-Ay n + — Az n 
ox Oy Oz 



(12) 
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From these equations we can easily infer that the components of dvi/dxj can be calculated from 



/ dv x \ 
dx 
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dv z 
~W 
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Av y2 

V Av y3 J 
fAv zl \ 

\Av z3 J 



(13) 



where A 1 is the inverse of the matrix 



A = 



Ax 2 Ay 2 Az 2 



(14) 



\Ax 3 Ay 3 Az 3 J 

Note that the four points of the interpolation tetrahedron are both necessary and sufficient to fix 
the value of each of the 9 velocity gradients. While in the linear regime only 6 of these quantities 
would be needed, by virtue of the absence of vorticity, all 9 quantities are necessary in the nonlinear 
regime. From the values of dvi/dxj we can then easily evaluate the value of the velocity divergence 
6, the shear cr^- and the vorticity uj (see eqn. 10) in each Delaunay tetrahedron. 

Subsequently we have to determine the corresponding volume averaged quantities. As de- 
scribed in the former section we accomplish this by top-hat filtering with a filter Wth that has a 
characteristic radius R. The problem has therefore been reduced to determining the average value 
of 9, dij or uj in a sphere of radius R, centered on some location xo- For example, for the volume 
averaged velocity divergence, 6>(x), we have 



0(xo) = 



J dx#(x)WV(x,x ) 



(15) 



As the value of 9 is constant within each cell of the space-filling Delaunay tessellation, the problem 
reduces to simply determining the intersection of the Delaunay tetrahedra with the filter sphere 
(appendix B), and using the intersection volume as the weighting value of 9 in integral of equation 
(15). In other words, 9 follows from, 

3 



0(xo) = 



47^ 



k 



(16) 
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where the sum is over all Delaunay tetrahedra k that intersect with the filter sphere, and Vk,R is 
the volume of the corresponding intersections, while 6k is the value of 6 in k. It is good to note 
that the geometrical issue of determining the intersection between a tetrahedron and a sphere is 
far from trivial. A description of our algorithm is provided in Appendix B, but we should note that 
more efficient evaluations are probably possible. Evidently, the discussion is exactly analogous for 
the volume-averaged values of the shear tensor components, of the vorticity, Cuij. 

As with the Voronoi method, the final end result of the operations described above consists of 
a field of top-hat averaged quantities at regular grid intervals. 

4. Computational Considerations 

In the former section we developed the theoretical groundwork for the use of the Voronoi and 
the Delaunay tessellations as optimal tools for determining the velocity divergence field, as well as 
the vorticity and shear fields, from the value of the velocity field at a discrete number of locations. 
The practical implementation of both methods consists of a two-phase approach. (1) In the first 
step, the algorithm calculates the Voronoi or Delaunay tessellation for the set of points at which the 
velocity is known/measured. From the Voronoi or Delaunay networks we can then define the value 
of the velocity divergence, vorticity and/or shear throughout the whole of the sample volume. (2) 
The subsequent second phase is the volume filtering of the resulting tessellation defined field. In 
essence, volume filtering consists of the determination of the average value of the velocity gradient 
field, produced by the tessellation algorithm, within the filter volume. Usually, the filtering averages 
are determined for a set of regular grid positions. 

4.1 The tessellation algorithm 

On the basis of the duality between Voronoi and Delaunay tessellations we use the same tes- 
sellation algorithm for both the Voronoi and the Delaunay method. It is the geometric Voronoi 
tessellation code that was developed by Van de Weygaert (1991b, 1994). For a very extensive and 
detailed description of the code we refer to both these references. Here we limit ourselves to a short 
overview of the fundamental ideas. 

The central operation of the geometric Voronoi algorithm is the calculation of all the Delaunay 
tetrahedra of the point distribution. Evidently, this simultaneously yields all Voronoi vertices, the 
centres of the corresponding circumscribing spheres. Finding the appropriate connections between 
the vertices, such that they define the edges, the polygonal Voronoi walls and finally the complete 
polyhedral Voronoi cells, is the major and most cumbersome part of the work. Our algorithm is 
an elaboration on the work of Tanemura et al. (1983), who gave a sketch of the basic algorithm 
for finding a Delaunay tetrahedron and a Voronoi polyhedron, in combination with four underlying 
geometric theorems. Through efficient administrative, data storage and searching procedures our 
algorithm manages to restrict itself to calculating every Delaunay tetrahedron and Voronoi wall only 
once during the construction of the complete tessellation. In a brute force construction procedure 
wherein each Voronoi cell is completely calculated without use of prior information this would 
be repeated three and one time respectively. This has the additional advantage that it allows to 
discard a point i from any further consideration once its Voronoi cell IL. has been calculated. 
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For the sake of clarity it is good to recapitulate some basic properties of Voronoi tessellations 
(also see Fig. 2). A nucleus i defines a complete polyhedral Voronoi cell. The boundary of the cell 
consists of a set of polygonal walls, each of which is shared with a nucleus j that is contiguous to i. 
A Voronoi wall is the part of space that is as close to i as to j, and closer to these two nuclei than to 
any of the other nuclei. On average there are ~ 15.54 walls per Voronoi cell if the generating point 
distribution is Poissonian, an average number that can vary from realization to realization (unlike 
the 2-D case, where the average number of edges is always 6). The boundary of a polygonal wall is 
formed by a number of edges, with a Voronoi vertex at the two tips of each edge. In other words, 
the structure of a wall is completely determined by listing the locations of the Voronoi vertices 
along its edge in the right geometrical order. Note that each edge is shared by three contiguous 
nuclei i, j and k, to which every point on the edge is equidistant while being closer to any of these 
three than to any of the other nuclei. Likewise, a Voronoi vertex is equidistant to its closest four 
nuclei in the generating point process. For a Poisson point process there are ~ 27.07 vertices per 
Voronoi cell, and some 5.228 per Voronoi wall. 

Our code assumes that the generating point processes, and consequently also the tessellations 
themselves, have periodic boundary conditions. This assumption, however, can be relaxed with 
some additional effort. Usually the points are distributed within a box of equally sized edges, 
although any rectangular parallelepiped is feasible. The first step of the algorithm is the storage 
of the set of generating nuclei in a multidimensional binary tree. This data structure, on average, 
stores points that are close together in space at nearby positions of a binary tree structure (see 
e.g. Bentley 1986). Nearby points in space can then be tracked efficiently by means of a recursive 
searching procedure. Building this tree is an N log N procedure, finding the nearest M points is 
a task that requires a limited effort proportional to M log N (see appendix F in Van de Weygaert 
1994). 

The algorithm then proceeds with the sequential computation of each Voronoi cell LTj. In turn, 
the construction of a cell IT consists of the determination of the structure of each of its Voronoi 
walls. A slight complication is of course that beforehand it is unknown with which and with how 
many nuclei i shares a wall. The identity of these contiguous nuclei can be revealed in different 
ways. Firstly, a subsample Si of M « 50 nearest nuclei is selected. As contiguous points are likely 
to be close in physical space, the subset Si is considered to be the primary sample of contiguity 
candidates. Initially any search for a new contiguous nucleus is restricted to Si. If % was not yet 
part of a Delaunay tetrahedron, the first contiguous nucleus is the nucleus %\ that is closest to i. 
Evidently, i\ belongs to Si. A subsequent second contiguous nucleus %2 is found by looking for the 
nucleus that together with i and i\ yields the triangle of minimal circumradius. Although almost 
always %i E Si, there are occasional exceptions, and an appropriate correction procedure for this 
has to be included. A subsequent third contiguous nucleus is the one that together with i, i\ 
and %2 defines a Delaunay tetrahedron, the first one of which i is one of the defining nuclei. As 
in the case of %2, the search for 13 is initially restricted to Si, but occasionally corrective action is 
needed and the search extended to a limited set of nuclei outside Si. This corrective procedure is 
efficiently performed via the multidimensional binary tree. On the other hand, i can have been one 
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of the four denning nuclei of one or more previously determined Delaunay tetrahedra. If so, we 
know that the other three nuclei of any of these tetrahedra are contiguous to i. In both situations 
an initial set Cj of contiguous nuclei is produced. Subsequently, we work our way through the list 
Cj, calculating one by one the corresponding walls. For each nucleus i a in Cj at least one Delaunay 
tetrahedron is known. These tetrahedra are the starting point of any further search. They are 
ordered in the appropriate geometrical order, and the gaps are filled in. Within the wall around 
{i,i a } a tetrahedron T>p a = {i, i a , ip, i 7 } should connect to a tetrahedron Vpj = {i, i a , ip, is}. If 
such a tetrahedron Vpj already exists its centre will be the Voronoi vertex connecting to the one 
of Vp~ r If not, we have to search for the nucleus is- First in the subset Si, if necessary followed 
by a corrective procedure. In this way new Delaunay tetrahedra are not computed at random, but 
always starting from previous knowledge of three of the four defining nuclei. It will automatically 
yield a new contiguous nucleus, and is is added to the list Cj. Continuing along the same direction, 
the polygonal wall gets progressively mapped by the connecting Voronoi vertices, a process which is 
completed once a tetrahedron connects to T>p^ on the side {i,i a ,i^}. Meanwhile newly determined 
Delaunay tetrahedra have produced new contiguous nuclei for the set Cj. After completion of 
the wall around {i,i a } the process proceeds with the next contiguous nucleus in the list Cj. The 
construction of the Voronoi polyhedron Ilj has been completed once the walls around all the nuclei 
in the constantly updated list Cj have been determined. 

In the construction procedure that we sketched in the previous paragraph each Delaunay tetra- 
hedron, and its corresponding Voronoi vertex, is computed only once. Once four nuclei i, i a , ip 
and z 7 are found to constitute a Delaunay tetrahedron, this information is passed on to all four 
of these nuclei. Likewise with the Voronoi wall shared by the nuclei i and j. During construction 
of the Voronoi cell Ilj the elaborate wall construction procedure sketched above is skipped when 
the wall with j has been computed earlier during the construction of Ilj. Instead, construction 
proceeds with the next nucleus in the contiguity list Cj for which not yet such a wall has been 
determined, naturally only after having notified the cell ilj of the existence of the wall it shares 
with j. Besides selection of the subset Si of M closest neighbours, one of the first steps in the com- 
putation of a Voronoi cell Ilj is therefore processing of the information on the walls and vertices 
that were obtained previously during the computation of other cells. This also implies that the 
nucleus i can be discarded from any further considerations once the cell Ilj has been determined, 
all of the Delaunay tetrahedra in which it partakes already having been calculated. This results in 
a progressive pruning of the multidimensional tree during the tessellation construction procedure. 

4.2 The filtering algorithms 

The tessellation algorithms lead to constant values of the velocity gradients within either the 
Voronoi walls, for the Voronoi method, or Delaunay tetrahedra, for the Delaunay method. Sub- 
sequently, this field is filtered (see Section 3.1 and 3.2). The operation of filtering is basically 
equivalent to determining a weighted average of the field within the filter volume. In the applica- 
tion of the Voronoi and the Delaunay method we restrict ourselves to determining top-hat filtered 
fields, which has the virtue that its filter value is a bounded sphere of radius R. 

Because in the Voronoi method the value of the velocity gradients differs from zero only in 
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the walls of the Voronoi tessellations, the filtering consists of determining the surface area of the 
intersections of the polygonal Voronoi walls with spheres of radius R. The technicalities of this 
procedure are treated in Appendix A. 

In the Delaunay method we end up with the more complicated situation of constant non-zero 
values of the velocity gradients in the Delaunay tetrahedra. The basic operation of top-hat filtering 
is therefore determining the volume of the intersections of tetrahedra with a sphere. This turns 
out to be a far from trivial geometrical problem. In Appendix B we describe a provisional recipe 
to accomplish this operation. However, it is quite likely that considerably more efficient ways are 
feasible, which would lead to a substantial gain in efficiency of the Delaunay method. 

4.3 Computational effort 

In practice, the amount of CPU time for the calculation of a Voronoi cell in the Voronoi 
tessellation construction procedure is almost independent of the number N of generating nuclei, so 
that the total tessellation construction time is almost linearly proportional to N. The construction 
of the multidimensional binary tree is a procedure that demands kN log N time. In addition, a time 
oc AM log N is required for the selection of each subset Si of M closest neighbours requires, so that 
a total time A(iVM) log A is spent on selecting all sets Si. It is the subsequent step for the actual 
construction of each polyhedral Voronoi cell IL, that dominates the CPU consumption. It is almost 
independent of N, due to the fact that it takes place in subsets Si that contain approximately the 
same number of M ~ 50 nuclei. Its total CPU requirement is therefore something like fiN, where 
H 2> A ~ k. 

To illustrate the above, we quote some CPU tests on an IBM 7011/25T Power PC. For a 
range of 100 to 50,000 generating nuclei the CPU time per cell is indeed constant, 0.010 sec. 
For a tessellation of 40,000 nuclei the tessellation computation time is « 407.28 sec, while the 
corresponding tree building time is « 28.03 sec. For comparison, in the case of a tessellation of 
1000 Voronoi cells this is 10.05 sec and 0.25 sec respectively. 

The construction of the Voronoi and Delaunay tessellation is considerably less time consuming 
then the ensuing Voronoi wall or Delaunay tetrahedron intersection procedures. The major share 
of the effort lies in the computation of the intersections between the spheres and either the Voronoi 
walls or Delaunay tetrahedra. The CPU time per intersection is basically constant. Here we quote 
some numbers for a top-hat filter of R = 15.0/j" 1 Mpc, where the filtered quantities were determined 
at 20 3 grid positions in a box with an edge size of 200. 0h~ l Mpc, containing a tessellation generated 
by A = 40,000 nuclei. The number of intersection evaluations is proportional to R 3 , as well as 
linearly proportional to N. With respect to the latter it is worthwhile to realize that a Voronoi 
tessellation generated by N nuclei has approximately 7. 768 A Voronoi walls and 6. 768 A Delaunay 
tetrahedra (i.e. ss 310, 720 walls and ss 270, 720 Delaunay tetrahedra for a Voronoi tessellation of 
40,000 cells). 

In the case of the Voronoi method we need on average 0.126 msec per intersection. The surface 
area of approximately 5,759,500 intersections had to calculated, which resulted in a total CPU time 
of 1454 sec. As expected, an intersection between a sphere and a Delaunay tetrahedron is more 
time demanding, taking « 1.32 msec, while the subsequent determination of the velocity gradient 
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Figure 6. Comparison of the estimated values of the velocity divergence 8 at the present epoch, top-hat filtered with a filter 
radius R = 15h~ 1 Mpc, at 20 3 grid positions. The simulation consists of 128 3 particles in a periodic box of 200/i~ x Mpc size, 
and concerns an Q = 1 Universe with CDM initial conditions (see text). The left-hand panel is a scatter plot of the value of 8 at 
each of the grid positions obtained with the two-step method against the value at the same position for the Delaunay method. 
The right-hand panel is a similar scatter plot, but then for the values obtained with the Voronoi method and the Delaunay 
method. 

Figure 7. Scatter plots of the value of the top-hat filtered velocity divergence 8 determined by the Voronoi method against the 
value determined by the Delaunay method, for 20 3 grid positions. The plot concerns the same 128 3 CDM N-body simulation 
as in figure 6. Each frame represents a different top-hat filter radius. Left-hand panel: R = 2.5k" 1 Mpc. Central panel: 
R = 5.0k- 1 Mpc. Right-hand panel: R = lO-O/i" 1 Mpc 

tensor (eqn. 13) takes on average another « 0.157 msec. Our code needed 12,882,264 intersection 
determinations, so that some 17,034 CPU seconds (~ 4.7 CPU hours) were spent on this part of 
the procedure, with an additional 1235 seconds for the velocity gradient tensor calculations. Hereby 
we should add the remark that due to the way the code selects the spheres and tetrahedra that are 
likely to intersect, some 39% of the intersections actually turned out to be empty. 

An additional important computational consideration concerns memory space. Here the De- 
launay method has a clear advantage over the Voronoi method. For the Voronoi method we need 
to store all the available information on the Voronoi tessellation. In this way we are able to recon- 
struct the Voronoi walls, in particular the links between and locations of the vertices that delineate 
those walls. As yet we still use the output of the general Voronoi code, although we intend to 
make a special-purpose version of the Voronoi code. The latter should reduce the memory alloca- 
tion considerably. For example, the structure of the complete Voronoi tessellation of 40, 000 cells 
is contained in 7 files of in total 105 Mbyte. By contrast, the Voronoi code was adapted for the 
Delaunay method so that only the information on the Delaunay tetrahedra is stored, comprising 2 
files of in total 14.3 Mbyte. 

In practice the limitations of memory space are considerably more restrictive than the required 
CPU time for the amount of data points that can be handled by the Voronoi and the Delaunay 
tessellation methods. For example, in the case of a data set of 128 3 particles, the Delaunay method 
would require a feasible 750 Mbyte of memory space. However, the more than 5.5 Gbyte needed by 
the Voronoi method probably prohibit this method of becoming applicable to data sets comprising 
more than « 100, 000 points. 

5. Application: Velocity Statistics of an N-body simulation 

In order to compare the new 'Voronoi tessellation method' (section 3.1), the new 'Delaunay 
tessellation method' (section 3.2) and the old 'two-step method', which we shortly described in 
section 2 (see e.g. Juszkiewicz et al. 1995, Bernardeau 1994b), we have applied the new methods to 
the result of an N-body simulation. This simulation, kindly provided by H. Couchman, was obtained 
with an adaptive P 3 M code (Couchman 1991) with CDM initial condition for Hq = 50Mpc/km/s, 
and periodic boundary conditions in a cubic box of 200/i~ 1 Mpc size. The simulation consists of 
128 3 particles and has been evolved until the epoch when the r.m.s. linear density fluctuations 
reaches unity in a spherical top-hat box of radius 8/i _1 Mpc. 

As it was not feasible for us to apply the Voronoi method for a particle sample of 128 3 points 
(see discussion at the end of section 4), we needed to reduce the number of particles used in the 
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analysis. To do so we take advantage of the fact that the filtered velocity field at the radius 
corresponding to the quasi linear regime is not sensitive to the details of the small scale velocity 
field in the very dense regions. Therefore, the statistical quantities we are interested in are not 
affected if we sample the data set in such a way that it does not reduce the number density of 
particles in the underdense regions while it does so in the very dense areas. To achieve this goal 
we have constructed an algorithm in which eight different grids are superposed, all with the same 
grid size but shifted by half of the grid-size in each of the possible directions. Subsequently, we 
perform a sequential operation on the set of simulation particles, whereby we consider the location 
of each particle, rejecting those whose locations in all of the shifted grids would be in a grid-cell 
that has already been occupied by a previously considered particle. For a grid size of lO/i^ 1 Mpc 
the final number of points is about 40000. It corresponds to a mean distance between points of 
about 6 /i _1 Mpc. We constructed the Voronoi and Delaunay tessellations of these particle samples, 
in order to derive the statistical properties of the velocity field in the way that was described in 
details in the previous section. In addition, we checked the robustness of the results by repeating 
the same analysis for a grid-size of 7.5/i _1 Mpc and double amount of points. 



Table 1. Statistical parameters for the density and the velocity field at different radii 
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1.50 ±0.12 


T-iVor. 


1.70 ± 0.05 


1.64 ± 0.07 


1.56 ±0.09 


1.48 ±0.12 


T 4 


5.35 


4.04 


3.25 


2.47 


y£>el. 


4.9 ±0.6 


4.6 ±0.6 


4.5 ± 1.5 


4.5 ±2 




4.6 ±0.4 


4.3 ±0.9 


4.2 ± 1.5 


4.3 ±2 



5.1 Results for the velocity divergence 

In order to illustrate the substantial improvements that are obtained by the Voronoi tessellation 
method and the Delaunay tessellation method, in comparison with older methods like the two-step 
one, and to compare the results to known theoretical predictions, we will in particularly focus our 
analysis on the velocity divergence. The same simulation was tentatively analyzed with the two- 
step method by Bernardeau (1994b, also see Juszkiewicz et al. 1995). In the scatter plot of figure 
6 we compare the estimates of 6 that we obtained with the Delaunay method on 20 3 locations 
on a regular grid, with those obtained by the two-step method (left-hand frame) and the Voronoi 
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Figure 8. The probability distribution function (PDF) of the top-hat filtered velocity divergence 6 at the present epoch, for 
the same CDM 128 3 particle N-body simulation as in figure 6 and 7, determined from the values on a 50 3 grid (see text). Top 
panel: filter radius Rth = 10/i -1 Mpc. Bottom panel: filter radius Rth = 15/i -1 Mpc. The solid lines represent theoretical 
predictions of the PDF for the measured values of ag (Bernardeau 1994b). The top panel one is given by eqn. (22), see 
Bcrnardeau (1994b), while the one in the lower panel has been obtained by numerical integration (see Bernardeau 1994b). 
Black squares: the Voronoi method. Black triangles: the Dclaunay method. Open squares: two-step grid method. 

method (right-hand frame), both on the same grid locations. The scatter plat clearly shows the 
good agreement between the Voronoi and the Delaunay tessellation method. Given the fact that 
both methods are quite different, this provides considerable confidence in them. On the other hand, 
the old two-step method yields very noisy estimates. Moreover, it tends to underestimate the value 
of by a factor of about 1.2. This effect is probably due to the fact that the effective smoothing 
radius is slightly larger because of the combination of the two smoothing procedures. 

In figure 7 we present similar scatter plots, for three different smoothing radii, to show the 
correlation between the divergences measured by the Voronoi and the Delaunay methods. The 
noise becomes very important for radii smaller than 5/i _1 Mpc, that is when the radius becomes 
smaller than the mean distance of the sampled particles. The results obtained for a smoothing 
length of 2.5/i~ 1 Mpc are obviously not reliable and show specific features associated with the 
methods that have been used. For example, at these small radii a large fraction of the smoothing 
spheres do not intersect any of walls of the Voronoi tessellation. The Voronoi method therefore 
yields values of zero for in these spheres. This situation is actually rather similar to the case of 
Poisson noise, with the velocity divergence only having a non-zero value at some discrete locations. 
For the Delaunay method the effect is less dramatic. However, the measured velocity divergence 
gets affected by the fact that information from larger scales, i.e. from the mean separation scale, 
leaks into the local velocity. This effect can also be traced in the other scatter plots. In the latter 
we note that the measured values of the velocity divergence tend to be smaller in the case of the 
Delaunay method than in the case of Voronoi method. Overall, however, at a scale of 10-15/i _1 Mpc 
we expect all measured quantities to be free of systematic errors, so that they can be analyzed with 
confidence. 

In Table 1 we summarize the resulting values of the moments of the velocity divergence obtained 
with the Delaunay method (with superscript DeL ) and the Voronoi method (with superscript Vor ), 
from an analysis in which the velocity divergence has been measured at 50 3 different grid points. The 
error bars have been obtained by dividing the simulation in four equal subsamples, and performing 
the same analysis for each of the subsamples. The results on the variance clearly show that the 
r.m.s. of the velocity divergence, ag, is lower than the r.m.s. density fluctuations, a§. This effect 
was also noticed by Juszkiewicz et al. (1995), and is an indication for the departure of the dynamics 
from the linear approximation, since in the linear regime ag and as are expected to be equal. 

However, our main interest concerns higher order moments. In the past years a lot of theoretical 
results on the statistical properties of the velocity divergence have been obtained. Bernardeau 
(1994a) derived the following results on the third and fourth order moment for cosmological models 
with Gaussian initial conditions, = 1 and A = 0. For the condition of a small ag, the following 
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expressions were found, 

^=T 3 + 0(aj(R)) (17) 

with 

26 

T 3 = — + 71 , (18) 

and 

{eA) -jf )2 =T 4 + 0(a$(R)) (19) 

with 

12088 338 7 2 2 
^ = ^ + -71 + 37^37,, (20) 

where the parameters 71 and 72 are the successive logarithmic derivatives of <r 2 (i?) with respect to 
R, 



dlogR 

dHogal{R) [ZL) 
72 dlog 2 R • 

For a CDM spectrum, Table 1 lists the values of the parameters n and 72 for different smoothing 
radii. The corresponding theoretical values of T3 and T4 can be determined from equations (18) 
and (20). The values for T3 measured from the simulation by both the Voronoi and the Delaunay 
method are found to be in remarkable agreement with these theoretical predictions in the case 
of all four different radii. While it was not possible to determine T4 as accurately in the case of 
the largest radius, it was found to agree reasonably well with the theoretical predictions for radii 
R^, 10/i -1 Mpc, i.e. at radii for which this quantity could be measured sufficiently accurately. This 
solves the issue raised by Lokas et al. (1995), who questioned the validity of the perturbation theory 
for the divergence of the velocity field. From our results we can conclude that the departure they 
observed in their work was due to the systematic errors introduced by the smoothing schemes they 
used. 

5.2 The Probability Distribution Function (PDF) of the velocity divergence 

Although it is interesting to study the individual moments, as they highlight different features 
of the total distribution, a more complete picture is obtained by looking at the global shape of the 
PDF of the velocity divergence. In figure 8 we present, for two different radii, the measured velocity 
divergence PDF. These measures are compared to the theoretical predictions of Bernardeau (1994b) 
obtained from the re-summation of the series of the cumulants (solid curves) . For n = — 1 there is 
a simple analytical expression that can be used for the PDF of 9 (given here for 17 = 1), 

P (9)d9 = 



([2k - lj/n 1 / 2 + [A - 1J/A 1 / 2 )- 3 / 2 



K 3/4 (2^)1/2^ 

with 



2Act 2 



< 22 > 



q2 off 

1 + — , and A = l-y. (23) 
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Figure 9. Comparison of the values of the top-hat averaged vorticity ui (eqn. 24) (left-hand frame) and shear a (eqn. 25) (right- 
hand frame) for the same CDM 128 3 particle N-body simulation as in figures 6-8, for a top-hat filter radius of Rth = 15/1" 1 Mpc 
(both u> and a in units of H). Both frames represent scatter plots of the values of these quantities at 20 3 grid positions as 
determined by the Voronoi method against those determined at the same position by the Delaunay method. 

This expression was used in the plot for R = lOft^Mpc. For R = 15ft, -1 Mpc, however, the solid 
curve was obtained by numerical integration of the inverse Laplace transform (similar to Eq. 18 of 
Bernardeau 1994b). 

The PDFs measured by both the Voronoi and the Delaunay method are clearly in remarkably 
good agreement with the theoretical predictions, down to probabilities of about 10~ 4 . On the 
other hand, previous methods, like the two-step method (open symbols in Fig. 8) produce PDFs 
that deviate substantially from the theoretical curves, producing spurious tails at both the low and 
high value end of the PDF. In particular noteworthy is the fact that the new tessellation methods 
manage to reproduce the expected sharp cutoff at the positive values of 9, which corresponds to 
voids. 



Table 2. The relative magnitude of the divergence, vorticity and shear. 



Radius/(/t _1 Mpc) 


7.5 


10. 


12.5 


15. 




0.63 ± 0.02 


0.53 ±0.02 


0.45 ±0.015 


0.38 ±0.015 




0.67 ±0.02 


0.55 ±0.02 


0.46 ± 0.015 


0.39 ±0.015 


^Del. 


0.183 ±0.003 


0.118 ±0.001 


0.081 ±0.001 


0.059 ±0.001 




0.243 ± 0.004 


0.147 ±0.002 


0.098 ±0.0015 


0.070 ± 0.001 


<7 Del - 


0.50 ±0.01 


0.41 ±0.01 


0.347 ± 0.009 


0.297 ±0.008 


ZfVor. 


0.55 ±0.01 


0.45 ±0.01 


0.366 ± 0.009 


0.310 ±0.008 



The accuracy of the perturbation theory calculations is therefore confirmed for the shape of 
the PDF of the velocity divergence. This is clearly of great importance, as the shape of the PDF 
can be potentially used to measure (Bernardeau et al. 1995). 

5.3 Vorticity and Shear 

When introducing the Voronoi and Delaunay method in section 3, we made the observation 
that both methods can actually be used to study the statistical properties of any quantity related 
to the velocity deformation tensor, such as the vorticity and the shear. In figure 9 we compare by 
means of scatter plots the results that have been found with our two methods for the norms of 
these quantities, uj (left-hand frame) and a (right-hand frame), 

- 2 = E-L (24) 
k 

<J 2 = J2 a v a ir ( 25 ) 

These plots should be compared with the similar plot of the velocity divergence in the right-hand 
frame of figure 6. Because the mean vorticity is pretty small its statistics is quite sensitive to noise. 
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Figure 10. Log-log plots of the probability distribution functions (PDFs) of the vorticity ui (in units of H, left-hand frame) 
and shear a (in units of H, right-hand frame) of the same 128 particle CDM N-body simulation as in Fig. 6-9, determined 
from the values on a 50 3 grid. Both frames concern the values of these quantities smoothed with the same filter as in figure 9, 
a top-hat filter with radius of Rth = 15/i — 1 Mpc. Black squares: Voronoi method. Black triangles: Delaunay method. 

Figure 11. Scatter plots of the value of the local density contrast 5 against the velocity divergence 8 (left-hand frame), vorticity 
u> (central frame) and shear a (right-hand frame), at 20 3 grid positions, for the same 128 3 particle CDM N-body simulations 
as in Figs. 6-10. All quantities are top-hat filtered with a top-hat radius Rth = 15/i _1 Mpc. The solid line in the left-hand 
panel indicates the prediction of the 5 — 9 relation by Bernardeau (1992). 

Even for a radius as large as 15/i _1 Mpc the resulting measured value of the mean vorticity can be 
significantly affected by systematics errors. Since it does not vanish in the linear regime such a 
sensitivity does not exist for the velocity shear. A summary of the expectation values of u and a 
is given in table 2, which also contains the values obtained for the r.m.s. of 6. 

We find that the amount of shear is slightly smaller than the amount of divergence. In fact, 
the value for a found with both methods is almost exactly consistent with the magnitude of the 
shear expected in the linear regime, 

jf 2 {<? 2 ) = I (0 2 ) • (26) 

As far as the amount of vorticity is concerned, we see that it shows the expected rapid decrease with 
scale. For R ^ 10/i _1 Mpc it seems to be fair to assume that the vorticity is small compared to the 
divergence. To give a crude idea of them, figure 10 shows the PDF of w and <r, for R = 15/t _1 Mpc, 
obtained with the two new methods. 

5.4 The local density-velocity relationship 

In addition to the analysis of the statistical properties of the velocity field, it is also possible to 
use the Voronoi and the Delaunay method to study the joint distributions of the density and the 
velocity field. Figure 11 displays scatter plots of the local density contrast 5 against respectively 
divergence (left-hand frame), vorticity (central frame) and shear (right-hand frame). The strong 
correlation between the density and the divergence is as expected, although there is a quite a large 
amount of scatter. For comparison, the solid line shows the prediction by Bernardeau (1992). From 
the scatter plot against the vorticity we can infer that the mean vorticity increases slightly with the 
density. A similar conclusion can be made for the shear. It confirms the idea that voids tend to be 
regular spherically expanding regions, whereas dense matter concentrations tend to have non-radial 
motion. 

6. Summary and Discussion 

The velocity field in the local Universe is an important and essential source of information on 
structure formation. In particular interesting for an understanding of the evolution and dynamics 
of the structures in the Universe are the various components of the gradient of the velocity field, the 
velocity divergence, shear and vorticity. One approach is to study the statistical properties, both 
moments as well as the full probability distribution function (PDF), of the divergence, shear and 
vorticity of the local smoothed velocity field. Considerable effort has been directed to obtaining 
analytical results for the statistics of the velocity divergence in the linear and quasi-linear regime 
in the case of structure formation scenarios based on Gaussian initial density and velocity fields. 
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To study more advanced stages of structure evolution we often have to resort to N-body sim- 
ulations, yielding discretely sampled velocity fields. The discrete nature of the velocity sampling 
complicates the determination of the statistical properties of the velocity field. In this paper we 
have introduced and developed two numerical methods, the Voronoi tessellation method and the 
Delaunay tessellation method, that yield reliable and accurate estimators of volume-averaged quan- 
tities in the case of discretely sampled velocity fields. The fact that they concern volume-averaged 
quantities is of crucial importance. Almost all analytical results concern volume- average quanti- 
ties while in essence all available numerical estimators only yield mass-averaged quantities. The 
latter considerably obscured the comparison between statistical results from analytical models and 
N-body studies, and even lead to false conclusions regarding e.g. the validity of perturbation theory. 

The availability of estimators of volume-averaged velocity statistics is important for several 
reasons. Firstly, it allows us to check independently whether the perturbation calculations that 
yield the quasi-linear results are indeed valid. Secondly, if so, we can apply the new estimators 
with confidence to highly nonlinear circumstances. And finally it may be feasible to apply them, 
in adapted form, to the available catalogues of measured galaxy peculiar velocities. 

The Voronoi tessellation method and the Delaunay tessellation method. Both methods are 
based on two important objects in stochastic geometry, the Voronoi and the Delaunay tessellation 
of a point set. A Voronoi tessellation of a set of nuclei is a space-filling network of polyhedral cells, 
each of which delimit the part of space that is closer to its nucleus than to any of the other nuclei. 
The Delaunay tessellation is also a space-filling network of mutual disjunct objects, tetrahedra in 
3-D. The four vertices of each Delaunay tetrahedron are nuclei from the point set, such that the 
corresponding circumscribing sphere does not have any of the other nuclei inside. The Voronoi and 
the Delaunay tessellation are closely related, and are dual in the sense that one can be obtained 
from the other. 

In a first evaluation of the two new methods we calculated, on a regular grid, the volume- 
averaged velocity divergence, shear and vorticity of an 128 3 particle N-body simulation of structure 
formation in an = 1 CDM Universe. Computer memory space limitations in the case of the 
Voronoi method forced us to sample some 40, 000 particles from the total sample of 128 3 particles. 
This sampling was performed such that the number density in underdense regions was not reduced. 
A comparison study between the Voronoi method, the Delaunay method, the conventional 'two- 
step' method and analytical theoretical predictions yields encouraging results for our new methods. 
The comparison study consists of comparison of scatter plots, third and fourth order moments as 
well as the global PDFs. The Voronoi and the Delaunay method show remarkable good agreement 
with each other, as well as with theoretical predictions. On the other hand, considerable differences 
with the conventional 'two-step' method were found. 

We may therefore conclude that the Voronoi and the Delaunay methods represent optimal 
estimators for determining the probability distribution function of volume-averaged velocity field 
quantities. As yet it is difficult to judge which of the two methods is the preferable one. The 
results produced by both methods agree very well for a considerable range of situations. However, 
the Delaunay method is clearly the better one for small filter radii, as it provides a reasonable 
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estimate of the velocity gradients throughout the whole of space. The Voronoi method, on the 
other hand, only does so in the Voronoi walls. Consequently, irrelevant and noisy filter averages 
are produced by the latter if the filter scale is smaller than the average wall distance because the 
small filter spheres frequently end up being empty. Another advantage of the Delaunay method is 
its approximately 8 times lower memory space requirement. However, this may be a mere practical 
issue, a more efficient Voronoi method implementation is certainly feasible. A clear disadvantage 
of the Delaunay method is the fact that it is almost 12 times more CPU time consuming than 
the Voronoi method. This is for a large part due to the inefficient calculation of the intersection 
between a tetrahedron and a sphere. We expect that better and faster prescriptions are possible, 
which may possibly lead to a five to tenfold acceleration of this algorithm. 

In forthcoming work we will apply the newly developed tool to a plethora of structure formation 
scenarios, based on both Gaussian as well as non-Gaussian initial conditions. The reliability of the 
results obtained with both the Voronoi and the Delaunay method allows us to study in how far 
the velocity field PDFs are sensitive discriminators that highlight physical differences between the 
scenarios. In particular, we are interested in the possibility of extracting the value of f2 from these 
velocity statistics. The availability of the reliable numerical estimators that we developed here is a 
crucial step in making this a practical possibility. We therefore also wish to develop our methods 
for the even less ideal circumstances of observational data. To see whether this is feasible, one of 
the first steps is to see in how far our methods are sensitive to substantial amounts of noise in the 
data. These issues will be addressed in a forthcoming study. Also, we intend to make the software 
that we developed for both the Voronoi and the Delaunay method publicly available once it has 
been made user-friendly. 

In addition, we feel that variations of the two numerical tools that we have introduced here 
can be applied to a choice of other applications in astrophysical situations. In many situations 
the value of a particular physical quantity is only known at a limited number of discrete points in 
space. The optimal adaptive nature of both the Voronoi and the Delaunay tessellations to the point 
distribution makes methods based on them promising estimators of the general run of quantities 
over the whole of the sample space. 

Acknowledgments We are very grateful to Dick Bond, Christophe Pichon and Simon White 
for useful and encouraging discussions and suggestions, and to Hugh Couchman for providing the 
results of the CDM N-body simulation. In addition, we would like to thank Bhuvnesh Jain for useful 
suggestions, and Jens Villumsen for encouraging remarks on Voronoi tessellations. In particular 
we are indebted to the referee, Edmund Bertschinger, for useful suggestions and comments. FB is 
grateful for the hospitality of the MPA, and RvdW for the hospitality of CE de Saclay, where most 
of the work presented in this paper was carried out. 

References 

Bentley, J.L., 1986, Programming Pearls, Addison- Wesley 
Bernardeau, F., 1992, ApJ, 390, L61 
Bernardeau, F., 1994a, ApJ, 433, 1 



25 



Bernardeau, F., 1994b, A&A, 291, 697 

Bernardeau, F., Juszkiewicz, R., Dekel, A., Bouchet, F.R., 1995, MNRAS, in press 
Bertschinger, E., Dekel, A., 1989, ApJ, 336, L5 

Bertschinger, E., Dekel, A., Faber, S.M., Dressier, A., Burstein, D., 1990, ApJ, 364, 370 
Burstein, D., Davies, R.L., Dressier, A., Faber, S.M., Stone, R.P.S., Lynden-Bell, D., Terlevich, 

R.J., Wegner, G.A., 1987, ApJS, 64, 601 
Coles, P., 1990, Nature, 346, 446 
Couchman, H.M.P., 1991, ApJ, 368, L23 
Dekel, A., 1994, ARAA, 32, 371 
Dekel, A., Rees, M., 1994, ApJ, 422, LI 

Delone, B.V., 1934, Bull. Acad. Sci. (VII) Classe Sci. Mat., 793 
Dirichlet, G.L., 1850, J. reine angew. Math., 40, 209 
Finney, J.L., 1979, J. Comput. Phys., 32, 137 
Gilbert, E.N., 1962, Ann. Math. Stat., 33, 958 

Goldwirth, D.S., Da Costa, L.N., Van de Weygaert, R., 1995, MNRAS, in press 
Icke, V., van de Weygaert, R., 1987, A&A, 184, 16 
Ikeuchi, S., Turner, E.L., 1991, MNRAS, 250, 519 

Juszkiewicz, R., Weinberg, D.H., Amsterdamski, P., Chodorowski, M., Bouchet, 1995, ApJ, 442, 
39 

Kiang, T., 1966, ZfA, 64, 433 

Lokas, E.L., Juszkiewicz, R., Weinberg, D.H., Bouchet, F.R., 1995, MNRAS, in press 

Matsuda, T., Shima, E., 1984, Prog. Theor. Phys., 71, 855 

Meyering, J.L., 1953, Philips Res. Rept., 8, 270 

Miles, R.E., 1970, Math. Biosc, 6, 85 

M0ller, J., 1989, Adv. Appl. Prob., 21, 37 

Nusser, A., Dekel, A., 1993, ApJ, 405, 437 

Peebles, P.J.E., 1980, The Large-Scale Structure of the Universe, Princeton Univ. Press 
Rubin, V.C., Thonnard, N., Ford, W.K., Roberts, M.S., 1976, AJ, 81, 719 
Smoot, G.F., Lubin, P.M., 1979, ApJ, 234, L83 

Stoyan, D., Kendall, W.S., Mecke, J., 1987, Stochastic Geometry and Its Applications, Akademie- 
Verlag, Berlin 

Strauss, M.A., Davis, M., Yahil, A., Huchra, J.P., 1990, ApJ, 361, 49 

SubbaRao, M.U., Szalay, A.S., 1992, ApJ, 391, 483 

Tanemura, M., Ogawa, T., Ogita, N., 1983, J. Comp. Phys., 51, 191 

Van de Weygaert, R., 1991a, MNRAS, 249, 159 

Van de Weygaert, R., 1991b, Ph.D. thesis, Leiden University 

Van de Weygaert, R., 1994, A&A, 283, 361 

Van de Weygaert, R., Icke, V., 1989, A&A, 213, 1 

Villumsen, J.V., 1995, MNRAS, submitted 

Voronoi, G., 1908, J. reine angew. Math., 134, 198 



26 



Figure 12. The intersection of a sphere and a polygon. The dashed disc is the intersection of the sphere with the plane 
supporting the polygon. The points I, J are the intersecting points of the polygon with the small circle. 

Figure 13. Sketch of the calculation of the volume Vaj 7 
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Appendix A: The intersection of a sphere and a polygon 

For the implementation of the Voronoi method we need to calculate the surface of the inter- 
section of a polygons and a sphere. 

Figure 12 sketches the geometrical situation we encounter. The problem naturally reduces to a 
planar problem in the plane of the polygon (bottom panel) where one has to calculate the surface 
of the intersection of the polygon with the disc. This figure corresponds to the generic case where 
the circle intersects the polygons in two different points I, J. The surface is then given by the sum 
of the surface of the polygon (I,A,B,C,D,J) and the surface of the segment defined by the points 
I, J, given by ir(9 — sm(9))R^/2, where is the radius of the disc. There are in fact other possible 
cases: the disc may be entirely contained in the polygon (the surface is then the one of the disc), 
may contain the whole polygon (the surface is the surface of the polygon) or may have more than 
two intersection points with the polygon. In the latter case the surface is obtained by a combination 
of polygons and segments. 

Appendix B: The intersection between a sphere and a tetrahedron 

For the implementation of the Delaunay method we need to calculate the volume of the in- 
tersection of a sphere and a tetrahedron. In the following we use (A,B,C,D) to indicate the four 
points that define the tetrahedron. The letters I, J, K and L represent any arbitrary order of these 
four points. 

The volume of the intersection is calculated by a sequence of complementary volume calcula- 
tions. To be specific, the intersection volume follows by taking the volume of the whole sphere as 
a start. From that volume we then extract the volumes cut out by each of the planes defined by 
three points (I,J,K). In total there are four of such planes, the possible permutations of (A,B,C,D). 
Subsequently, we have to correct by adding each of the volumes contained in the six spherical 
segments defined by two planes (I,J,K) and (I,J,L). Finally, we should subtract the volumes of the 
four cones defined by the three planes containing either I, J, K, or L. In short, 

^intersection — 

Sphere ~ Vp l,J,K + J2 Va i,J ~ S Vs i ' 

perm. perm. perm. 

where the summations are made over the possible permutations for I,J,K,L, and 
• ^sphere is the volume of the sphere; 
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Figure 14. Sketch of the calculation of the volume V$j 

• Vpi ik ls * ne volume of the sphere segment delineated by the plane (I,J,K) on the side opposite 
to the point L. 

• Va, j is the volume of intersection of the sphere segments carved out by the planes (I,J,K), 
opposite to the point L and (I,J,L), opposite to the point L; 

• Vsj is the volume of the intersection of the sphere segments defined by the planes (I,K,L), 
opposite to J, (I,J,L), opposite to K and (I,J,K), opposite to L. 

The calculation of Vp is quite straightforward. It is given by the expression, 

V p = itR 3 (2/3- x-x 3 /3), (B2) 

where x is the distance of the center of the sphere to the plane in units of the radius, R, or by its 
complementary part in the sphere. 

The calculation of Va t , is intrinsically more complicated. The geometrical problem is illus- 
trated in Figure 13 in the plane orthogonal to the planes (I,J,K) and (I,J,L). The edge (I, J) is 
indicated by the point /. The distance x is the distance of the centre of the sphere to this line 
(expressed in units of the radius). The volume to be calculated, Va, is indicated by the doubly 
shaded area. It depends on x and the angle 9, and for xtan# < 1 is given by 

V A = 

- \—dx 2 sin 9 cos9 
3 1 



+ x (3 + (3 — x 2 ) tan 2 9) sin 9 cos 2 9 arctan ^- 



+ arctan 



— arctan 



d 2 tan9 



d(l — x tan i 
x + d 2 tan 9 



x cos# N 

) 

(B3) 



d 



d(l + xtan0) 
+ x it (3 sin 3 9 + 3 cos 2 9 sin 6* - x 2 sin 3 (9)/2} 
and by the same expression plus 7r/3 otherwise. In this expression, the parameter d is defined by 

d= (1 -x 2 ) 1/2 . (B4) 

This expression is strictly valid only when < 9 < ir. Otherwise eqn. (B2) and (B3) have to be 
combined to get the proper answer. 

The calculation of the volumes Vs is even more cumbersome for its practical implementation, 
but is given by a combination of equations (B2) and (B3). The typical geometrical situation is 
presented in Figure 14. By default we assume the tip I to be located within the sphere. If this 
is not the case the situation is simpler and can be reduced quite straightforwardly to previously 
discussed situations and equations. 

The volume Vsj can also be calculated by the determination of a sequence of complementary 
volumes, 

V Sl = Vtetra. " £ Vp, + £ V Aj , K „ (B5) 
perm. perm. 
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where Vtetra. is the volume of the tetrahedron defined by I and the intersection points J',K',L' of 
the three half- lines [I, J), [I,K) and [I,L) with the spheres. The volume Vpi is the one of the fraction 
of the sphere above the plane (J',K',L') (obtained with eqn. B2), and Vaj, k , are the three volumes 
of the fractions of the sphere that are bounded by (J',K',L') and respectively (I,J',K'), (I,J',L'), 
(I,K',L'), and that do not contain I nor respectively L', K', J'. These volumes are given by a proper 
use of expression (B3). 
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